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ABSTRACT 


A  fourier  transform  method  is  used  to  invert  the  integral 
equation  resulting  from  holographic  interferometry  of  three- 
dimensional  weakly  refractive  phase  objects.   The  inverted 
equation  is  applied  to  the  determination  of  aerodynamic 
density  fields.   The  effects  of  various  numerical  techniques 
are  extensively  evaluated  in  an  attempt  to  minimize  computa- 
tional effort.   Numerical  filtering  techniques  are  introduced 
in  order  to  handle  the  discontinuities  associated  with 
supersonic  aerodynamic  phenomenon. 
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I.   INTRODUCTION 

Holographic  interferometry  offers  a  means  of  determining 
three-dimensional  density  fields  created  by  aerodynamic 
phenomenon.   Pulsed  lasers,  with  good  spatial  and  temporal 
coherence  with  double-exposed  holograms  have  eliminated  the 
need  for  high  quality  optical  components.   The  short  pulse 
of  the  giant-pulse  lasers  allows  the  aerodynamicist  to 
"freeze"  transient  phenomenon  in  the  interferometric  views 
around  the  test  section. 

Interferometric  data  is  obtained  by  superimposing  a 
hologram  of  the  test  section  without  flow  and  a  hologram  of 
the  test  section  with  flow.   The  superposition  is  accomplished 
by  double  exposure  of  the  same  holographic  plate.   When  the 
double-exposed  hologram  is  viewed,  phase  differences  in  the 
two  views  result  in  light  and  dark  regions.   By  counting 
the  number  of  light  and  dark  regions  in  going  from  an  un- 
disturbed region  to  the  point  of  interest,  one  obtains  the 
fringe  number  for  that  point.   The  phase  difference  is 
determined  by  the  total  difference  in  optical  path  length, 
resulting  from  the  aerodynamic  phenomenon.   The  relative 
refractive  index  [f(x,y)]  is  a  function  of  the  density 
[p(x,y)]. 

f(x,y)  =  -^-  Cp(^y)  -1)  Eq.  (1) 

PSA       Poo 


where  p   is  the  density  of  the  gas  at  which  g  is  measured. 
s 


To  determine  the  refractive  index  from  the  fringe  data 
available  around  the  test  section  g(p»8)j  one  must  invert 
the  integral  equation. 

'  L2 
g(/>,0)  =   /  f(x,y)  d» 

By  utilizing  a  two-dimensional  fourier  transform  method 
[Ref.  *1],  it  has  been  shown  that  the  refractive  index  at 
a  point  (X  ,Y  )  in  the  test  section  can  be  expressed  as 


2 


f(x^,Yj  -  -iy  .  /   /  g(po+p,e>+s(P0-P,e)-2g(po,e)d 

n 

"2 


j      j      v         ' u    * _u u.  p  d  0 

2tt    tt     n  2 

0  p 


The  fourier  transform  method  requires  less  computational 
effort  than  the  previously  utilized  methods  of 


1.  Fourier  expansion  of  the  unknown  function  f(x,y) 
in  appropriate  orthogonal  polynomials.   [Ref.  1] 

2.  Representing  the  unknown  functions  by  discrete 
cubes  and  solving  the  resulting  set  of  simultaneous 
linear  equations. 


The  computerized  method  developed  yields  results  which 
compare  favorably  with  the  fourier  expansion  method  of  Matulka 
[Ref.  1].   The  computation  time  was  reduced  by  a  factor 
greater  than  twenty.   Errors  resulting  from  the  fourier 
transform  method  were  greater  than  those  associated  with  the 
fourier  expansion  method  of  Reference  1.   The  reduced  compu- 
tation time  more  than  offsets  the  slightly  larger  errors 
associated  with  the  fourier  transform  method. 
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II.   HOLOGRAPHIC  INTERFEROMETRY 

The  ability  of  a  hologram  to  record  both  amplitude  and 
phase  variations  has  greatly  extended  the  capabilities  of 
interferometric  analysis.   When  a  hologram  storing  more  than 
one  wave  pattern  is  illuminated  with  coherent  light,  the 
reconstructed  waves  interfere  with  each  other,  producing 
fringe  patterns  which  yield  quantitative  information 
concerning  the  disturbance  causing  the  interference.   Since 
the  two  wavefronts  were  constructed  utilizing  the  same 
optical  paths,  the  optical  components  are  automatically 
matched  and  only  the  changes  in  the  subject  create  the 
visible  fringe  patterns. 

The  effective  reconstruction  of  the  wavefront  by  holo- 
graphy requires  that  the  coherent  light  of  a  laser  be  sepa- 
rated into  two  beams  (Pig.  1),  the  object  or  scene  beam 
which  illuminates  the  subject  (modulated),  and  the  reference 
beam  which  is  reflected  from  an  optically  flat  mirror  (un- 
modulated).  The  path-lengths  of  the  two  beams  must  be 
approximately  equal  due  to  the  finite  coherency  length  of 
available  lasers  (10  cm  -  1  m) . 

In  order  to  simplify  the  development  of  the  reconstruction 
process,  assume  the  reference  beam  is  collimated,  although 
this  is  not  necessary  in  practice.   Since  the  exposure  time 
of  the  film  represents  a  time  average  of  the  intensity,  the 
temporal  portion  of  the  wavefronts,  eJ w  may  be  neglected 


and  only  spatial  variations  in  the  plane  of  the  photographic 
plate  considered.   Let  U  (x,y)  denote  the  complex  amplitude 
of  the  reference  beam.   Since  this  beam  is  plane 


Uo(x,y)  =  ao  e1(ox  +  6y) 

where  a  and  3  are  the  spatial  frequencies  of  the  reference 
beam  in  the  x,y  plane  of  the  holographic  plate  (cycles/mm) 

Similarly  let  U(x,y)  denote  the  complex  amplitude  of 
the  scene  beam 


U 


(x,y)  =  a(x,y)  e1^*^ 


The  intensity  I(x,y)  thus  recorded  by  the  photographic 
film  is  given  by  the  expression 

I(x,y)  =  |U+U  |2  =  a2+a  2faa  e*[*<*.y)-<«-f»y] 

1    o  ■        o    o 

+aa  e~1^^X5y^"ax~ey^ 
o 

2    2 
=  a  +a  +2aa  cos[<J)(x,y)-ax-3y] 

When  the  exposed  hologram  is  developed,  the  trans- 
mittance  will  be  proportional  to  the  intensity  I(x,y). 
Illumination  of  the  hologram  by  a  plane  beam  U  (x,y)  similar 
to  the  reference  beam  yields  a  transmitted  beam  U.(x,y) 
given  by 
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Ut(x,y)  =  Ur(x,y)[l-kl(x,y)] 


=  [l-k(a2+ao2)]Ur(x,y)-ao2kU(x,y)-ka(x,y)U*2 


where  k  is  the  constant  of  proportionality  of  the  photographic 
emulsion. 

The  illuminating  beam  produces  a  direct  beam  and  two 

first-order  diffracted  beams  on  either  side  of  the  direct 

2 
beam.   The  term  ka  U(x,y)  represents  the  diffracted  beam 

that  produces  the  virtual  image.   The  last  term  produces 

the  real  image  with  inverted  depth.   The  negative  sign  for 

the  two  diffracted  beams  simply  represents  a  180  degree 

phase  shift  to  which  the  eye  is  insensitive. 

If  we  consider  only  the  virtual  image  formed  by  double 

exposure  to  two  complex  wavefronts,  the  transmitted  virtual 

image  is  given  as 

Ut(x,y)  =  -k[ao2(a1(x,y)ei*l(x'y)+a2(x,y)ei*2(x^))] 

Hence  illumination  of  the  double-exposed  hologram  not  only 
effects  the  simultaneous  reconstruction  of  the  two  waves 
but  causes  them  to  interfere. 

Although  it  is  not  necessary  for  the  original  reference 
beam  to  be  collimated,  if  the  re-illumination  beam  is 
different  from  the  reference  beam,  the  reconstructed 
diffraction  beams  are  distorted  accordingly.   Since  a 
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collimated  beam  is  the  easiest  to  duplicate,  a  collimated 
reference  beam  should  be  used  when  possible. 

For  a  more  detailed  and  complete  analysis  of  holography, 
Reference  2  is  recommended. 
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III.   FRINGE  INTERPRETATION 

By  utilizing  different  methods  of  superimposing  the 
two  wavefronts,  various  fringe  patterns  may  be  observed. 

A.  LIVE  FRINGE 

Live  fringe  patterns  are  generated  when  the  hologram  is 
exposed,  developed,  replaced  exactly  in  its  original  posi- 
tion, and  re-illuminated.   The  stored  wave front  in  the 
hologram  interferes  with  the  "live"  wavefront  from  the  sub- 
ject.  This  technique  requires  precise  repositioning  of  the 
hologram  plate  and  readjustment  of  the  scene  beam's  intensity. 

B.  FROZEN  FRINGE 

Frozen  fringe  patterns  are  produced  by  the  double  expo- 
sure of  the  hologram  prior  to  development.   The  double- 
exposed  hologram  need  not  be  restored  to  its  original  position 
in  order  to  view  the  fringe  patterns.   Illumination  by  a 
beam  similar  to  the  original  reference  beam  is  quite  adequate. 

C.  INFINITE  FRINGE 

Infinite  fringe  patterns  are  seen  if  the  two  wavefronts 
viewed  through  the  hologram  utilized  the  exact  same  optical 
arrangement  for  construction.   The  visible  fringe  patterns 
are  focused  at  infinity. 

D.  FINITE  FRINGE 

Finite  fringe  patterns  occur  when  finite  displacements 
are  introduced  in  the  optical  arrangement  between  exposures 
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of  the  hologram  plate.   The  displacements  usually  used  are 
small  rotations  of  one  of  the  mirrors  in  either  the  reference 
or  scene  beams,  translation  of  the  diffusely  illuminating 
screen,  or  translation  of  the  diffusely  illuminating  screen, 
or  translation  of  the  hologram  plate  [Ref.  3].   If  there  is 
no  disturbance  of  the  subject  between  exposures,  these  small 
translations  (d  =  rG  ~  .003  inches)  result  in  parallel  inter- 
ference lines  which  are  perpendicular  to  the  direction  of 
translation.   The  parallel  lines  form  a  convenient  scale  for 
determining  fringe  orders.   The  fringe  pattern  will  be 
focused  somewhere  within  the  phase  object  [Ref.  3]. 

In  conventional  holography,  three-dimensionality  is 
accomplished  by  recording  the  wavefront  created  by  the 
diffusely  reflected  light  from  the  subject  [Fig.  2].   Each 
point  on  the  subject  is  a  point  light  source  which  interferes 
with  all  other  point  light  sources  on  the  subject.   This 
same  three-dimensionality  may  be  retained  in  holograms  of 
transparent  objects  by  placing  a  diffuser  plate  in  the 
scene  beam  prior  to  the  test  section  [Fig.  1].   The  locali- 
zation of  fringe  patterns  is  greatly  enhanced  by  using 
these  "light-field"  interferograms . 

The  utilization  of  frozen,  finite  fringe  patterns  on 
lightfield  interferograms  yields  good  results  when  dealing 
with  aerodynamic  phenomenon  associated  with  wind  tunnel 
experiments . 
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IV.   FOURIER  TRANSFORMS 

With  appropriate  restrictions  [Ref.  2],  optical  systems 
may  be  considered  linear,  space  invariant  systems.   As 
linear  space  invariant  systems,  they  may  be  analyzed  using 
Fourier-transform  techniques.   Temporally  modulated  waves 
may  be  analyzed  in  either  of  two  domains;  a  temporal  domain 
or  a  temporal  frequency  domain.   Analogously,  spatially 
modulated  waves  may  be  analyzed  in  either  a  spatial  domain 
or  a  spatial  frequency  domain. 

In  the  spatial  domain,  the  light  complex  amplitude  a(x,y) 
is  expressed  as  a  function  of  the  (x,y)  spatial  coordinates. 
The  same  complex  amplitude  can  be  expressed  in  terms  of 
orthogonal  spatial  frequencies  a  and  g  (rads/mm)  as  well. 
Thus  the  light's  complex  amplitude  distribution  a(x,y) 
in  the  spatial  domain  may  be  expressed  as  another  function 
A(a,8)  in  the  spatial  frequency  domain.   The  two  functions 
are  related  through  the  two-dimensional  transform  pair 

A(a,B)  -  -~z   f°     /a(x,y)  ei(ax  +  6y)  dx  dy 
i/2tt  -< 


.00     _00 


and 


(x,y)  =  -pz  f*     /°A(a,S)  e"i(ax  +  By)  da  d3 

^TT  -oo    _oo 
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To  every  function  in  the  spatial  domain,  there  usually 
is  a  corresponding  function  in  the  spatial  frequency  domain. 
Similarly  for  every  operation  in  the  spatial  domain,  there 
is  a  corresponding  operation  in  the  spatial  frequency  domain 
The  advantage  of  Fourier  analysis  is  that  the  operation  in 
the  spatial  frequency  domain  may  be  much  easier  than  the 
corresponding  operation  in  the  spatial  domain.   The  multi- 
plier (1/  /2"tt)  is  used  since  most  Fourier  transform  tables 
list  the  transform  pairs  such  that  the  total  multiple  is 
1 

p-TT. 
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V.   THE  INTEGRAL  INVERSION 

A.   THE  INTEGRAL  EQUATION 

The  information  available  from  the  hologram  [Fig.  3]  is 
in  the  form  of  fringe  numbers.   The  relative  phase  delay 
experienced  by  straight  line  rays  passing  through  the  phase 
object  determines  the  fringe  number  at  some  point  of  the 
hologram. 

If  n(x,y,z)  is  the  index  of  refraction  with  the  phase 
object  present,  and  n  (x,y,z)  is  the  index  of  refraction 
without  the  phase  object,  the  relative  phase  in  wavelengths 
is 

2  n(x,y,z)  -  n  (x,y,z)-, 
g(p,0,z)  =    /   [ , 2 J  dv       (ia) 

Ll 

where  (Lp  -  L, )  is  the  total  distance  of  travel  of  each 
parallel  ray  through  the  phase  object.   From  Figure  3  the 
relationship  of  the  coordinates  (p,v)  to  (x,y)  for  a  given 
ray  direction  (6)  is 

x  =  v  sine  +  p  cose  (2A) 

y  =  -v  cose  +  p  sine  (2B) 

The  refractive  index  of  gas  as  a  function  of  density  is 
given  by  the  first  two  terms  of  the  Taylor  expansion  (Liepmann 
and  Rushko,1957) . 
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n  =  1  +  3  p/ps 


where  p   is  the  density  of  the  gas  at  which  3  is  measured 
s 

For  air  at  zero  degrees  centigrade  and  760  mm  Hg 


p   =  .0012929  gm/cc 
3   =  .00291   for  deep  red  light 
Hence : 


p  A  n(x,y,z)  -  n  (x,y,z)-. 
p(x,y,z)  =  pQ(x,y,z)  +  —  [ ^ 

(3) 


The  wavelength  (X)  is  maintained  because  of  the  form  of 
Eq.  (1A) 


n(x,y3z)  -  n  (x,y,z) 
Let    f(x,y,z)  =  5 2 


By  choosing  the  (x,y)  plane  parallel  to  the  rays,  the 
z  coordinate  becomes  a  parameter  and  Eq.  (1)  can  be  written 

L2 
g.(p,e)  =   /  f  (x,y)  dv  (IB) 

Ll 

The  parameter  z  will  be  dropped  from  further  discussion. 

Previous  solutions  to  the  three-dimensional  refractive 
index  problem  utilized  Fourier  series  techniques  [Ref .  1] 
and  required  excessive  calculations.   The  method  proposed  by 
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Rowley  [4]s  as  modified  by  Junginger  and  Haeringen  [5], 
offers  good  results  for  a  minimum  of  calculations. 

B.   THE  FOURIER  INTEGRAL  SOLUTION 

The  limits  of  Equation  1A  can  be  extended  to  infinity 
because  f(x,y)  is  zero  outside  the  phase  object. 

00 

g(p,e)  =   /  f(x,y)  dv  (1C) 

—  00 

The  function  f(x,y)  may  be  described  by  the  Fourier  integral 
f(x,y)  -  —   /"   /°F(a,e)  e"i(ax  +  3y)  da  dg  (4A) 

x/27T  -< 


.00    — oo 


Substituting  equations  (2A),  (2B),  and  (4A)  into  (1C)  and 
integrating  first  on  v  gives 


B(p.e)  -  -p=   /"  /"p(a.B)  e-1P(aC08et6Bine)da  dg 

V27T  -00   -°° 


^givCasine-gcose)  dv  (5) 


The  Dirac-Delta  function  has  the  following  transform. 


*(t0-t)  =  ±       /°°eia)(to-t)  da)  (6) 


—  00 


Hence   the   last   integral   of  equation    (5)    is    simply 


,      iv(asin6-3cos9)  ,  0    , ,       ,    .    „         nS 

f  e  M  ydv   =    2TrS(asin6-Bcose) 


^(atane-B) 


I  cose 
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Because  of  the  properties  of  the  delta  function,  if  we  now 
integrate  over  J}  >   Equation  (5)  reduces  to 


iP: 


a 


g(p,6)  =  v^ff   /  F(a,B=atan6)  e'   cos6  ugse|   (7) 


Equation  (7)  has  the  Fourier  inverse 

a 
F(a,B=atan9)  =   L   -/o        f   g(p,G)e   |cos01  dp   (8) 

Although  Equation  (8)  and  (*JA)  may  be  used  to  determine 
the  refractive  index,  there  are  numerical  difficulties 
associated  with  the  infinite  integrals  of  Equation  (4A). 
The  problem  can  be  further  simplified  as  suggested  by  Junginger 
and  Haeringer  [53- 

For  a  given  point  (x  ,y  ),  substitution  of  Equation  (8) 
into  Equation  (4A)  gives 

n        co     oo     oo  -i- 2°L_   ax  _£y   ) 

f(x  ,y  )  =  -\       f        f        f   g(p,e)e(  lcosel    °   °  dpdadg 

O    O       I,  c- 

L\y:  _0O    —00     —00 

(4B) 


Introducing  polar  coordinates  in  the  Fourier  frequency  domain 


0  =  arctan  g/a 


.  _  (    2    2*1/2  _   a 


|  cose 
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From  Figure  (3) 


p   =  x   cos8  +  y   sin< 
oo        ''o 


Equation  (*JB)  becomes 


2tt  a  .      a 


tiff       oo  oo  /IP"! 5T-P     o\ 

rf  .     \  1  r  r  r        r        a\      (  COS0  OCOS0) ,    .    ,,,,, 

f(x    ,y    )    =   — 7T     J         J         J    g(p,9)e         '  '  dpkdKdi 


4tt2    0         0       - 


Separating  the   integral   over  theta  yields 

— -     OO       00  . 

f(xojyo)  =  "^2    /2  ;    ;  g(P,0)  eik(p-po}dp  kdk  d6 

2   u 

2      00       00 

+  -i-  /   /   /  g(p,6)  eik(p+po)  dp  kdk  di 
J   0   - 

Utilizing  the  fact  that  g(-p,8)  equals  gCpjO+Tr),  the  last 
term  may  be  written  as 

7T_ 

2   °°   °° 
..+  -i   /  /   /  g(p,6)  e"ik(p~po)  dp  kdk  d9 

-IT    0    -oo 


Recombining  the   two   terms   yields 

TT_ 

2  oo  CO 

f(x    ,y    )    =  — *       f      f        f  g(p,6)cosk(p-p    )    dp  kdk  di 

27T         _7T_        0         -oo  ° 

2  (4C) 


Performing  the  p  integration  by  parts  yields 
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n 

2     oo      CO     .  . 

f(xo,yo)  =  — ig-   /  /   /  ^  p?   sink(p-po)dp  dk  de 
2?r   _£  0   -oo    p 

2  (9) 


since  g(±°o,0)  =  0. 

The  integrations  over  k  and  p  can  nov;  be  interchanged  and 
the  k  integration  can  be  evaluated.   Equation  (9)  becomes 

H 

2ir   1   -00  o 

"2 

where  £?  stands  for  the  principal  value  of  the  following 

integral.   Using  the  definitions  of  the  principal  value  and 

taking 

ag(p,e)  _  aCg<P.B)-g(p0,e)] 
dp  dp 

we    find   integrating  by  parts   over   p 

2L 

2 


f(xo'yo)  =  ~~2      J      f  tg(p0+P,e)-g(po,e)]^- de 

(10A) 


2.TC  £ 

^2 


Pc 

00  K 


which  can  be  rewritten  as 


2L 

2    oo 


f(x  ,y  )  =  ~±?      f    f   [g(P  +p,e)+g(P  -p,e)-2g(P  ej]%de 

2ir   n   0  P 

2  (10B) 


A  Taylor  expansion  about  p  =  0  shows  that  the  integrand  in 
(10B)  is  well-behaved. 

s 
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Equation  (10B)  demonstrates  that  the  controlling 
parameters  for  determining  the  refractive  index  at  a  point 
are  the  relative  phase  delays  near  the  point. 
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VI.   NUMERICAL  METHOD 

A.   PROGRAM  PARAMETERS 

The  integration  (from  -  tt/2  to  tt/2)  over  theta  yields 
the  same  results  as  integrating  over  phi  from  0  to  7r[Fig.  3]. 
The  angle,  phi  equals  zero,  corresponds  to  the  direction  of 
the  parallel  rays  or  viewing  angle.   The  integration  over 
phi  is  used  in  order  to  make  the  computer  program  coordinate 
consistent  with  the  laboratory  coordinates . 

The  FORTRAN  IV  program  developed  assumes  that  all  holograms 
have  the  same  number  of  equally-spaced  data  points  and  that 
the  total  number  of  data  points  on  each  plate  is  odd.   A 
square  grid  is  established  over  the  phase  object  with 
spacing  determined  by  spacing  of  the  data  points  on  the 
holograms  at  0  and  90  degrees.   The  views  are  assumed  to  be 
equally  spaced  and  to  cover  180  degrees  of  view  for  asymmetric 
phase  objects.   Symmetry  reduces  the  necessary  field  of 
view  until  only  one  view  is  required  for  the  axi-symmetric 
case. 

To  visualize  how  symmetric  fields  may  be  reconstructed 
from  reduced  fields  of  view,  consider  a  phase  object  of 
unity  refractive  index  as  shown  in  Figure  14.   The  field 
is  symmetric  about  the  phi  equals  zero  axis  (x  axis).   The 
fringe  information  in  the  second  quadrant  is  exactly  the 
same  as  that  in  the  first  quadrant.   The  numerical  method 
developed  generates  additional  data  as  necessary,  depending 
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en  the  degree  of  symmetry.   The  180  degree  view  is  in  all 
cases  generated  by  the  computer  program  from  the  zero 
degree  viev;  by  copying  the  fringe  data  in  reverse  order. 
For  a  description  of  input  parameters  see  the  computer 
program. 

B.   THE  INFINITE  p  INTEGRATION 

Ideally,  g(p,$)  would  be  zero  for  p  greater  than  RMAX 
(the  radius  of  the  inversion  domain).   Actual  data  from  a 
wind  tunnel  will  in  some  cases  have  a  non-zero  fringe  shift 
at  RMAX.   The  method  devised  to  handle  the  case  of  non-zero 
fringe  shifts  at  RMAX  expands  the  fringe  data  using  a  cubic 
polynomial  which  matches  the  slope  at  RMAX  and  is  zero  at 
2  RMAX.   This  method  is  used  in  all  cases  since  if  g(p,$) 
is  zero  at  RMAX,  the  coefficients  for  the  cubic  polynomial 
are  zero.   For  the  general  case,  g(p,$)  is  zero  for   p 
greater  than  2  RMAX. 

The  numerator  of  the  integrand  of  Equation  (10B) 
[g(P0+P>4>)+g(p0-Pj$)-2g(po,<$)]  is  not  a  function  of  p  for 
p  >  3  RMAX  provided  that  |p  |  <  RMAX.   The  infinite  p 
integration  is  divided  into  two  regions 


1    ,*  r3RMAX  g(po+p^-)+g(po-p,?)-2g(po,j,)dpdj) 


f<VV  -  77  '   *  -2 

2tt   0   0  D 


_1_       *        "  g(Po+P,0)+g(Po-P,»)-2g(po><»)dp  d^ 
2tt2  0   3RMAX  p2 
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The  integrand  in  the  second  term  reduces  to 


g(P0+P,$)  +  g(p5-P,4>)-2g(po,<M     2g(po,c()) 

2  =      2 

P  P 


because  g(p  +  p,<J>)  +  g(p  -p,<f>)  is  zero  for  p  >_  3  RMAX.  Since 
g(p  ,<{))  is  not  a  function  of  p,  the  integration  over  p  may 
be  performed  and  the  second  term  becomes 


1      7T  2g(PQ,«J>) 

/     O   PMAV ^ 


2?  o    3  ™AX 


Recombining  the  two  terms,  one  has 


-,  tt        3RMAX  g(p   +p,<J>)+g(p   -p,<f>)-2g(p    ,<|>) 

f(Vo}  =  -  -h  f  c  ;        — 2 — ~ 2 — dp 

00  2v      0        0  p 


2g(PoScf)) 
+    3   RMAX        ]    d*  (10C) 


Figure  12  shows  a  typical  variation  of  the  argument  of 
Equation  (10B)  versus  p.   The  p  integration  from  zero  to 
3  RMAX  is  performed  using  Cote's  sixth  order  formulas  [Ref.  10] 

For  the  points  (x  ,y  )  corresponding  to  the  established 
grid  intersections  (Fig.  3)5  the  values  of  p   correspond  to 
available  data  points  only  for  the  views  where  phi  equals  0 
or  90  degrees.   For  all  other  views,  the  value  of  p  will 
lie  somewhere  between  two  data  points.   The  procedure  used 
calculates  the  value  of  the  inside  integral 
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3RMAX  g(pL+p,*k)+g(pL-p,<j>k)-2g(pL,cj>k) 

G(pL,<f>k)  =   /      " 2 "  dp 

0  p 

2g(pL,»k) 

3  RMAX 


at  each  data  point  on  each  hologram.   Hence  for  each 
g(pT  >40  there  corresponds  a  function  G(pT,4>,  )  which  is  the 
inside  integral  of  Equation  10B  evaluated  at  that  data  point. 
For  a  given  view  <f>  and  a  given  point  (x  ,y  ),  the  correspon- 
ding value  of  p  will  be  such  that  [p,  <  p   <  pT,, ].   A  cubic 
&  Ko  KL  —  Ko  —  KL+1 

polynomial  in  p  is  fitted  to  the  points  G(p,  ,,$,),  G(p, <}>,), 

G(PL+1>4)k) »  and  G(PL+2'*k^  and  the  value  of  G^p0'*k^  is 
obtained  by  evaluating  the  cubic  polynomial  at  p  . 

If  NP  is  the  number  of  data  points  for  each  view  and 
K  is  the  number  of  views,  the  inside  integral  must  be 
evaluated  (NP-K)  times.   For  each  of  the  (NP)   grid  inter- 
sections there  are  K  interpolations  for  G(p  j<t>v)  and  an 
integration  over  phi. 


C.   THE  INTEGRATION  OVER  PHI 

The  number  of  views  obtained  in  actual  experiments  are 
usually  quite  limited.   These  limited  views  are  the  deter- 
mining factor  in  the  accuracy  of  the  inversion  process. 
It  is  desirable  to  achieve  a  fairly  accurate  method  of  inte- 
gration over  phi  that  does  not  place  undue  restrictions 
on  the  experimenter. 
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The  integration  on  phi  proceeds  from  0  to  tt;  however, 
the  value  of  the  inside  integration  is  known  outside  this 
range.   Let  G(4>)  be  the  value  of  the  inside  integration 
and  A<J>  equal  the  spacing  of  the  views.   Then 

G(0   -    Acf>)    =    G(l80    -   A<J>) 
G(l80   +    A(f>)    =    G(A<Jj) 

By  using  overlapping  parabolas  through  each  succeeding 
data  point  with  the  first  parabola  centered  at  phi  =  0, 
the  next  parabola  centered  at  phi  -  A<J>,  etc.,  the  following 
coefficients  result 

f(x  7  j  -  .  *%  [2^1  +  o(A4)  ♦  ...  +  o(i8o-A#)4aiya)] 

which  is  the  trapezoidal  rule  for  integration. 

Hence  the  integration  over  phi  proceeds  quite  rapidly 
and  no  restrictions  are  placed  on  the  experimenter  as  to 
the  number  of  views  that  must  be  supplied,  except  that  A<f> 
must  not  be  too  large. 
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VII.   PROGRAM  TESTING 

A.   TEST  INVERSIONS 

The  computer  program  has  the  capability  of  generating 
the  fringe  shift  information  g(p,<|>)  for  an  arbitrary 
refractive  index  which  is  supplied  by  the  user.   This 
feature  may  be  used  to  test  the  ability  of  the  program  to 
invert  various  anticipated  density  fields.   The  following 
test  refractive  index  fields  were  used  in  MODE  1  (see 
computer  program) . 

1 .   Axisymmetric  Gaussian 

The  relative  refractive  index  supplied  was 


f(x,y)  =  e(-x2-y2) 


with  a  2  cm  inversion  radius.   [Figure  *J]  ( 

2.   Axisymmetric  Shock 

The  function  supplied  was  similar  to  that  which 
would  be  encountered  by  a  supersonic  cone  at  zero  angle  of 
attack. 

f(x,y)  =  -  .5  Vx2+y2  +2   R  <  1.5 

=0  R  >  1.5  [Figure  51 
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3.   Axi symmetric  Shock  With  Restricted  Data 

The  function  was  as  in  (2),  except  with  restricted 
data  which  doesn't  have  a  sufficiently  large  inversion 
radius. 

f(x,y)  =  -  .5  Vx  +y2  +  2       R  <  1.5 

=  .25  1.5  <  R  <  3.5 

=0  R  >  3.5 

[Figures  6  &  7] 

k .   Asymmetric  Gaussian 

The  relative  refractive  index  was  given  by 


2 
f(x,y)  =  (1  -  y2)  e"x       R  <  1.2 


=0  R  >  1.2 

[Figure  8] 

B.   INVERSIONS  OF  ACTUAL  DATA 

Reference  1  is  concerned  with  the  analytical  and  experi- 
mental analysis  of  an  axisymmetric  free-jet.   Further,  the 
free-jet  was  tilted  to  provide  a  test  field  which  was 
asymmetric  in  the  plane  of  the  solution.   Since  the  input 
data  files  are  available  in  the  reference,  the  same  data 
was  used  to  test  the  fourier  transform  method.   Three  sets 
of  data  were  used  (see  input  data  files). 


30 


1.  An  axisymmetric  section  at  z  =  . 05  cm   [Figure  9] 

2.  An  axisymmetric  section  at  z  =  . 5   cm    [Figure  10] 

3.  An  asymmetric  section  at    z  =  .5   cm    [Figure  11] 
The  axisymmetric  data  consisted  of  101  data  points  at 

various  z  locations.   Asymmetric  data  was  available  at 
only  one  position  (z  =  .5)  and  consisted  of  nine  views  with 
ten  degree  spacing  of  the  views,  and  with  60  data  points 
per  view.   The  views  covered  only  90  degrees  since  there 
exists  one  plane  of  symmetry  for  the  tilted  free-jet. 
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VIII.   DISCUSSION  OF  RESULTS 

A.   TEST  INVERSIONS 

MODE  1  of  the  computer  program  artificially  generates 
the  g  array  by  performing  the  line  integral  of  Equation  (IB) 
over  the  specified  relative  refractive  indexv [f (x,y) ] .   In 
order  to  reduce  quantizing  effects  of  the  numerical  integra- 
tion, it  was  necessary  to  divide  the  line  integral  into  5*np 
divisions  where  np  is  the  number  of  data  points  from  each 
view.   There  was  little  reduction  in  the  quantizing  errors 
when  higher  order  numerical  integration  techniques  were 
used;   therefore,  the  trapezoidal  rule  was  implemented. 

MODE  2  generates  the  g  array  from  the  values  of  the 
refractive  index  which  is  supplied  at  the  grid  points. 
Although  this  allows  the  use  of  test  density  fields  which 
cannot  be  specified  analytically,  the  quantization  errors 
in  the  g  array  strongly  influence  the  inversion  process. 
The  primary  reason  for  this  section  of  the  program  is  to 
regenerate  the  g  array  for  MODE  4  operation  (see  computer 
program) . 

For  the  axisymmetric  case  the  integration  over  4>  is  in 
five  degree  increments  utilizing  trapezoidal  integration 
as  previously  discussed.   For  the  asymmetric  case,  the 
views  are  specified  by  the  user. 

1 .   Axisymmetric  Gaussian  Inversion 

The  axisymmetric  Gaussian  function  tests  the  ability 
of  the  program  to  invert  smooth,  continuous,  axisymmetric 
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density  fields  [Fig.  4].   Using  101  data  points,  the  error 
in  the  inversion  was  everywhere  less  than  three  per  cent 
of  the  maximum  value  of  the  specified  function.   The  largest 
errors  occurred  near  the  boundaries  of  the  inversion  region. 
2.   Axisymmetric  Step  Function 

The  axisymmetric  step  functions  [Fig.  5  &  6]  test 
the  ability  of  the  program  to  invert  functions  similar  to 
those  associated  with  supersonic  flow  around  axisymmetric 
shapes.   Although  the  g  array  is  continuous,  it  has  a  dis- 
continuity in  slope  at  the  step  [Fig.  13].   In  the  continuous 
case,  the  integrand  of  Equation  (IOC)  would  become  infinite. 
Because  of  the  discrete  nature  of  the  numerical  method,  the 
integration  over  p  yields  a  finite  value  even  though  it  may 
be  quite  large  if  a  large  number  of  data  points  are  used 
[Fig.  13].   The  severe  properties  of  the  function 


G(p  ,4>)  =   /  [g(po+P,4>)+g(po-p><J>)-2g(po,<j>)]  £| 
0  p 


create  oscillations  and  over-shoots  in  the  inversion.   It 
was  necessary  to  introduce  an  averaging  process  for  the 
integration  over  $ 

p  +kAp 
G(P„,*)  =  ^-  I  S(p)  dp 


o-"    ^P 


o  p  -kAp 
o 


where  S(p)  is  the  third-order  polynomial  through  the  four 
closest  values  of  the  inside  integral  and  k  is  specified 
by  the  user  (0  <_  k  <_  1).   This  averaging  process  could  be 
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considered  a  low  pass  filter  which  rejects  the  higher  fourier 
frequencies.   Although  the  process  limits  overshoots  and 
oscillations,  it  also  slows  the  response  to  the  step  change 
in  refractive  index.   The  amount  of  filtering  is  optional 
to  the  user  depending  on  the  value  of  k  chosen.   There  is 
no  filtering  if  k  =  0  and  maximum  available  filtering 
occurs  at  k  =  1. 

With  the  exception  of  points  near  the  discontinuity, 
the  maximum  errors  occur  at  points  external  to  the  step 
where  the  maximum  error  was  six  per  cent  of  the  maximum 
value  of  the  refractive  index.   Interior  to  the  shock,  the 
errors  in  the  inversion  were  on  the  order  of  two  per  cent. 
The  reason  for  the  larger  errors  outside  the  step  is  readily- 
apparent  from  Figure  13.   For  a  point  (x  ,y  )  outside  the 
step  function,  the  phi  integration  proceeds  across  the  severe 
portions  of  the  function  G($).   If  the  point  (x  ,y  )  lies 
within  the  step  function,  the  integration  over  phi  never 
crosses  the  severe  portion  of  the  curve  [Fig.  15]. 

.  The  value  of  the  inversion  tends  toward  the  average 
at  the  discontinuous  points  in  the  function  supplied.   The 
amount  of  filtering  used  affects  the  rise  distance  and 
overshoot  of  the  inversion  at  the  shock  as  would  be  expected 
of  any  low-pass  filter. 

3.   Restricted  Radial  Coverage 

Figure  7  demonstrates  the  feasibility  of  the  smoothing 
process  discussed  in  Section  VI-A.   The  fringe  shift  g(p,<j>) 


34 


was  computed  from  the  specified  function  using  a  radius  of 
lJ  cm.   Only  the  data  for  a  radius  of  2  cm  was  used  as  input 
for  MODE  3  (g  array  supplied).   Except  for  the  region  near 
the  boundary  of  the  inversion  domain,  the  resulting  inver- 
sion showed  no  degradation  from  the  results  of  the  inversion 
using  the  entire  g  array. 

^ •   Asymmetric  Test  Case 

The  asymmetric  test  function  was  chosen  such  that 
it  contained  discontinuous  regions  as  well  as  regions  of 
negative  relative,  refractive  index.   Six  views  were  generated 
over  the  required  90  degrees  of  view  with  101  data  points 
per  view.   The  inversion  was  accurate  to  within  5  per  cent 
over  the  entire  domain,  except  near  the  discontinuity.   The 
inversion  showed  good  response  to  the  discontinuity  [Fig.  8]. 
The  0  degree  cross  section  shows  that  the  program  continuously 
underestimates  the  refractive  index.   A  test  conducted 
using  the  axisymmetric  function  of  Figure  H   as  an  asymmetric 
function  with  the  same  number  of  views  verified  that  this  is 
the  case.   Since  the  only  difference  was  the  number  of 
views  available  (11  versus  37)  for  the  integration  over 
phi,  the  phi  integration  method  is  believed  responsible. 

B.   INVERSION  OF  ACTUAL  DATA 
1.   Axisymmetric  Free- jet 

The  input  data  available  from  Reference  1  contained 
100  data  points  per  hologram.   Since  it  was  necessary  to 
have  an  odd  number  of  data  points,  an  additional  data  point 
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was  added  at  p  =  0  by  averaging  the  points  on  either  side. 
The  results  of  the  inversion  process  were  within  3  percent 
of  the  results  obtained  by  Matulka  [Fig.  9  &  10].   There 
was  evidence  of  a  continuous  under-estimation  of  the  density 
field  by  the  program.   The  effect  of  adding  one  additional 
point  was  evident  near  the  origin  where  the  maximum  differ- 
ence in  the  two  solutions  occurred.   The  present  program 
also  has  a  slower  response  at  the  boundaries  due  to  the 
filtering  used. 

2.   Asymmetric  Free- jet 

An  additional  data  point  was  added  to  the  available 
data  at  p  =  0  to  give  a  total  of  61  data  points  on  each  of 
the  nine  views.   The  addition  of  the  data  point  is  believed 
to  be  the  primary  cause  of  the  low  value  for  the  density  at 
the  origin  [Fig.  11].   The  procedure  of  adding  one  data 
point  at  the  origin  by  an  averaging  technique  affects  not 
only  the  value  at  the  origin,  but  also  points  near  the  origin 
since  it  is  the  relative  phase  delays  that  are  being  measured 
The  reason  for  the  low  value  of  density  at  points  near  the 
origin  is  not  caused  by  the  fact  that  the  interpolated  value 
for  the  fringe  shift  at  the  origin  is  grossly  in  error. 
The  error  is  introduced  by  the  fact  that  the  data  from  a 
60  point  spacing  is  read  into  a  61  point  spacing. 

To  visualize  how  this  occurs,  consider  the  1.5  cm 
radius  of  inversion  used.   The  grid  spacing  for  60  data 
points  would  be  3.0/59.   The  grid  spacing  for  6l  data  points 
would  be  3.0/60.   The  two  grids  correspond  exactly  at  ±RMAX 
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and  deviate  by  half  a  grid  spacing  at  the  center.   This 
distortion  of  the  resulting  inversion  can  be  visualized  as 
a  stretching  of  the  density  field  near  the  origin  which  is 
gradually  removed  by  compression  of  the  density  field  as 
one  moves  toward  the  boundaries  of  the  density  field. 
There  was  a  ten  percent  difference  in  the  two 
solutions  at  the  origin,  and  approximately  five  percent 
difference  at  the  boundaries  of  the  jet.   Again  the  filtering 
process  slowed  the  response  of  the  inversion  in  regions 
with  large  gradients. 
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IX.   CONCLUSIONS  AMD  RECOMMENDATIONS 

The  Fourier  transform  technique  represents  considerable 
savings  in  computational  effort.   With  appropriate  restric- 
tions, it  may  be  used  successfully  in  the  analysis  of 
discontinuous  density  fields.   The  simplicity  of  Equation 
(10B)  makes  the  controlling  parameters  readily  apparent. 
The  integrand  of  Equation  (10B)  is  related  to  the  rate  of 
change  of  the  slope  of  the  g  array.   It  has  been  demonstrated 
that  filtering  techniques  can  be  implemented  that  will 
handle  cases  involving  g  arrays  with  discontinuous  slopes. 
The  computer  program  written  has  numerous  options,  is  easy 
to  use,  and  requires  little  computational  time. 

The  ability  of  the  program  to  handle  discontinuous  g 
arrays,  such  as  might  be  encountered  in  studying  complex 
flow  patterns,  has  not  been  fully  documented.   It  is  expected 
that  the  filtering  method  used  will  handle  discontinuous 
g  arrays. 

Additional  investigation  of  numerical  methods  to  deal 
with  the  problems  of  discontinuities  is  warranted.   The 
importance  of  numerous  views  around  the  test  section  must 
not  be  underestimated.   Due  to  the  restrictive  nature  of 
wind  tunnels,  it  is  often  difficult  to  obtain  the  necessary 
number  of  views.   Any  method  which  would  enable  the  experi- 
menter to  reduce  the  number  of  views  or  decrease  the  necessary 
field  of  view  would  be  invaluable. 
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